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1. 


SUMMARY 


The  major  breakthroughs  obtained  in  the  development  of  the  Composite  Computer  Aided 

Design  (CAD)  tools  (1)  for  mixed  technology  integration  in  this  effort  include  the 

developments  of: 

•  FASTSTOKES,  which  is  up  to  3  orders  of  magnitude  faster  than  conventional 
numerical  methods  of  determining  drag  forces  on  geometrically  complicated  3- 
dimensional  structures. 

•  NeKxar  code  for  analysis  of  microflows.  NsKxar  code  includes  slip  models  for 
smooth  as  well  as  rough  walls. 

•  Force  coupled  method  for  bioparticulate  flows  as  well  as  techniques  and  codes  for 
stochastic  modeling  of  microflows. 

•  Molecular  codes  based  on  Direct  Simulation  Monte  Carlo  (DSMC)  and  molecular 
dynamics  techniques  for  micro  and  nanoscale  analysis  of  gas  and  liquid  flows. 

•  Finite  cloud  and  boundary  cloud  meshless  methods  for  efficient  interior  and  exterior 
analysis  of  microfluidics  and  MEMS. 

•  Multiscale  methods  combining  meshless  and  continuum  based  microfluidic 
simulation  tools  with  DSMC.  Investigated  fundamental  physical  phenomena  in 
microfilters  using  multiscale  simulation  tools. 

•  New  reduced-order  models  based  on  weighted  Karhunen-Loeve  decomposition 
technique  for  fast  dynamic  analysis  of  microfluidics. 

•  A  framework  for  coupling  SPICE3  with  fluidic  solvers  NsKxar  and  EOFLOW. 
Implemented  efficient  time  stepping  schemes  for  coupled  circuit  and  microfluidic 
device  simulation. 

Applications  of  this  software  have  resulted  in  the  efficient  computational  prototyping  of 

mixed-technology  microfluidic  components  and  systems  in: 

•  Microfabrication  and  packaging  technologies  for  microfluidic  systems  based  on 
silicon,  glass,  elastomer,  biodegradable,  and  liquid  crystal  polymer  materials. 

•  Demonstrated  a  micromachined  elastomer  micropump  for  remote  wireless  operation 
and  a  magnetic  micro  bar  mixer  for  mixing  fluids  in  channels  and  reaction  chambers. 

•  Design  rules  for  the  optimization  of  90°  and  180°  constant  radius  turns  in 
electrophoretic  systems.  A  novel  experimental  technique,  bleached  fluorescence 
visualization,  was  implemented  to  verify  these  predictions. 

•  Optimized  turn  geometry,  shown  to  reduce  the  dispersion  variance  to  less  than  1% 
for  a  simple  turn,  was  proposed. 

•  Design  of  several  complex  fluidic  devices  with  the  developed  design  software. 
Several  fabrication  cycles  have  been  eliminated  because  of  the  use  of  these  design 
tools. 

•  Initial  transition  of  design  tools  developed  as  part  of  this  project  to  commercial 
MEMS  CAD  companies. 
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2. 


INTRODUCTION 


The  objective  of  this  effort  titled  “Efficient  Computational  Prototyping  of  Mixed- 
Technology  Microfluidic  Components  and  Systems”  included  developing  fast  and 
efficient  computer-aided  design  tools  for  microfluidic  components  and  systems. 
Innovation  in  numerical  methods  and  algorithms  is  critical  as  we  work  towards  the 
development  of  advanced  CAD  tools  for  mixed-technology  microfluidic  devices  and 
systems. 

To  accomplish  this,  algorithms  and  techniques  for  coupled-domain  simulation  of  mixed- 
technology  MicroElectro  Mechanical  Systems  have  been  developed.  In  particular, 
algorithms  and  techniques  for  coupled  circuit/device  modeling  of  microfluidic  systems 
have  been  developed  and  tested.  Test  structures  and  devices  were  fabricated  to  verify  our 
models,  and  measurements  were  made  on  the  fabricated  devices  to  enable  computational 
prototyping  of  MEMS.  This  technology  was  then  transferred  to  the  commercial, 
university  and  government  researchers. 
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3. 


MICROFLUIDIC  SIMULATION  TOOLS 


Three  approaches  were  used  to  develop  microfluidic  simulation  tools  for  microfluidic 
energy  domains.  FASTSTOKES  focuses  on  fast  exterior  analysis  of  microfluidic  energy 
domains.  Also,  highly  accurate  spectral  methods  for  analysis  were  also  developed. 
Thirdly,  new  computational  methods  for  radically  simpler  analysis  of  microfluidic  energy 
domain  were  explored.  Summaries  of  this  work  are  presented  here,  and  additional  details 
can  be  found  in  references  1-18. 

3.1  FASTSTOKES 

FASTSTOKES  is  a  precorrected-Fast  Fourier  Transform  (FFT)  accelerated  fluid  analysis 
program  for  the  analysis  of  the  Stokes  equation  that  is  as  much  as  three  orders  of 
magnitude  faster  to  run  than  existing  finite-element  and  boundary-element  methods  when 
used  to  analyze  complicated  micromachined  structures.  FASTSTOKES  is  the  only  fluid 
analysis  program  capable  of  analyzing  entire  complicated  three-dimensional  designs. 
Developing  a  robust  numerical  engine  and  extending  the  analysis  capabilities  led  to  a 
large  number  of  algorithmic  and  physical  developments. 

3.1.1  Basic  Solver  Technology 

Several  new  algorithms  and  a  new  software  platform  were  developed  and  demonstrated 
for  the  precorrected-FFT  solver.  The  algorithmic  innovations  included  a  new  polynomial 
projection/interpolation  algorithm  and  a  new  approach  to  vector  convolution  that 
accelerated  the  code  by  nearly  a  factor  of  three. 

One  of  the  improvements  made  to  FASTSTOKES  was  to  reduce  the  number  of 
convolutions  required  to  compute  grid  velocities,  given  grid  forces.  In  the  classical 
approach  (Figure  1),  the  grid  velocities  were  computed  from  grid  forces  using  9 
convolutions  (a  total  of  18  Fast  Fourier  Transforms  (FFT’s)).  The  improved 
FASTSTOKES  approach  (Figure  2)  required  only  6  FFT’s. 
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Figure  1.  Classical  approach  to  computing  grid  velocities 
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Figure  2.  FASTSTOKES  approach  to  computing  grid  velocities 


The  first  Greens  function  independent  fast  solver  was  developed  and  demonstrated, 
which  makes  use  of  sparse  matrix  data  structures  to  make  it  easy  to  adapt  the  code  to 
other  applications.  The  code  has  been  used  for  applications  as  diverse  as  electromagnetic 
interconnect  analysis  and  offshore  structure  analysis  (FASTWAMIT),  as  well  as  for 
Stokes  flow  analysis. 

Also,  the  first  approach  to  computing  integrals  over  curved  panels  was  demonstrated  that 
is  both  accurate  and  also  independent  of  panel  separation  distance.  The  technique  used  a 
novel  implicit  mapping  scheme  to  map  the  curved  surface  to  a  flat  one,  a  recursive 
algorithm  for  integrating  singularities  over  flat  surfaces,  and  a  specialized  quadrature 
algorithm. 


3.1.2  Innovations  in  FASTSTOKES  Flow  Analysis 

The  integral  equations  for  Stokes  flow  have  a  well-known  singularity  for  single  body 
problems.  Because  the  codes  were  too  time  consuming  and  computer  intensive  to  use  on 
multiple  body  problems,  the  multidimensional  singularity  for  multiple  bodies  had  not 
been  previously  examined.  An  approach  to  eliminating  the  singularity  and  computing  the 
correct  drag  force  distributions  by  using  a  combination  of  a  singular-subspace 
orthogonalized  Krylov  subspace  algorithm  and  an  auxiliary  integral  pressure  boundary 
condition  was  developed. 


3.1.3  Application  of  FASTSTOKES  to  Comb  Drive  Analysis 

The  FASTSTOKES  Solver  was  used  to  compute  the  drag  on  comb  drives  (Figure  3),  and 
the  results  were  correlated  with  the  measured  data  from  Professor  Dennis  Freeman's 
(MIT)  visualization  system  for  MEMS.  The  detailed  drag  force  analysis  can  only  be 
made  by  using  FASTSTOKES,  as  the  drag  can  not  be  directly  visualized,  and  showed 
that  only  half  the  drag  force  on  a  comb  drive  is  due  to  comb- substrate  interactions.  This 
overturned  many  commonly  accepted  analyses  in  the  literature. 
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Figures.  A  comb  resonator  (dimensions  in  microns) 


3.1.4  Extensions  of  FASTSTOKES  to  the  Unsteady  and  Compressible  Case 

The  FASTSTOKES  code  was  extended  to  analyzing  both  unsteady  and  mildly 
compressible  Stokes  flow.  Several  new  techniques  were  developed  for  computing 
integrals  of  the  frequency-dependent  Greens  function  for  the  unsteady  case.  A  new 
linearization  approach  was  developed  that  can  analyze  compression  without  volume 
discretization  for  the  compressible  case.  The  latter  result  made  it  possible  to  analyze  the 
vertical  motion  in  micromachined  devices,  which  introduces  compression,  nearly  as  fast 
as  the  analysis  of  lateral  motion,  where  there  is  no  compression. 


3.2  SPECTRAL  METHODS 

Very  efficient  spectral/hp  methods  for  analysis  of  Stokes  and  Navier-Stokes  equations 
were  developed.  The  spectral/hp  methods  incorporate  both  multi-domain  spectral 
methods  and  high-order  finite  element  methods.  NsKxar,  a  spectral/hp  element  code, 
was  developed  for  microflow  analysis.  Specific  accomplishments  in  this  area  included 
the  development  of: 

•  Triangular/Tetrahedral  Spectral  Elements 

•  Slip  models  when  rough  walls  are  present  in  microfluidic  devices. 

•  A  force  coupled  method  for  bioparticulate  (2)  flows. 

•  Techniques  and  codes  for  stochastic  modeling  of  microflows. 

•  Molecular  codes  for  analysis  of  microflows. 
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3.2.1  N8Kxar 

NsKxar  is  a  general  purpose  Computational  Fluid  Dynamics  (CFD)  code  for  simulating 
incompressible,  compressible  and  plasma  flows  in  unsteady  three-dimensional 
geometries.  The  code  uses  meshes  similar  to  standard  finite  element  and  finite  volume 
meshes,  consisting  of  structured  or  unstructured  grids  or  a  combination  of  both.  The 
formulation  is  also  similar  to  those  methods  corresponding  to  Galerkin  and  discontinuous 
Galerkin  (3)  projections  for  the  incompressible  and  compressible  Navier-Stokes 
equations. 
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Figure  4.  Hierarchy  of  the  NsKxar  code  for  transient  nonlinear  flow  analysis 


In  Figure  4  is  presented  the  hierarchy  of  the  NsKxar  code  for  transient  nonlinear  flow 
analysis  (4,5).  The  “2d”  and  “3d”  refers  to  two-dimensions  and  three-dimensions, 
respectively.  “2.5d”  refers  to  a  three-dimensional  capability  but  with  one  of  the 
geometric  directions  being  homogeneous  in  geometry.  “ALE”  refers  to  moving 
computational  domains  required  in  dynamic  flow-structure  interactions.  Gaseous 
microflows  can  be  simulated  by  either  the  compressible  or  incompressible  version 
depending  upon  the  pressure/density  variations. 

Sets  of  numerical  simulations  of  dispersed  two-phase  flow  at  low  to  moderate  Reynolds 
numbers  for  both  unbounded  domains  and  ducts  with  complex  geometry  have  been 
performed.  The  code  NsKxar  has  been  employed,  with  the  3D  fluid  motion  fully 
resolved  by  the  spectral/hp  element  method  and  the  coupled  momentum  and  energy 
between  the  particle  and  fluid  phases  solved  by  finite-valued  force  coupling  method. 
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This  included  both  force  monopole  and  dipole  terms  together  with  inter-particle  and 
particle-wall  collision  barriers. 


Software  tools  for  post-processing  and  visualization  of  data  were  implemented  in 
MATLAB.  This  allowed  the  verification  of  the  results  for  solid,  fluid  and  circuit 
elements  of  Microsystems.  As  an  example  of  an  analysis  using  NsKxar,  flow  in  a  two- 
dimensional  microfluidics  channel  has  been  modeled.  As  shown  in  Figure  5,  the  two 
ends  of  the  channel  are  10  microns  long  and  1  micron  wide.  The  middle  channel  is  10 
microns  long  and  10  microns  wide.  The  fluid  in  the  middle  channel  obeyed  continuum 
theories,  and  deviations  from  classical  continuum  theories  was  observed  in  the  two  thin 
channels. 


1  a 


^  mem 

o 


Figure  5.  Two-dimensional  microflow  modeled  by  NsKxar 


3.2.2  Low-Dimensional  Vorticity-Based  Models 

On  the  physics  side,  out  objective  was  to  construct  low-dimensional,  vorticity-based 
dynamical  models  for  prototype  flow-structures  interaction  problems,  based  on  the 
simulation  databases  and  supplementary  experimental  information.  In  addition,  a  more 
fundamental  goal  was  to  provide  insight  to  a  selected  class  of  these  flows  and  explain 
such  generic  phenomena  as:  Large-scale  instabilities  of  non-stationary  turbulent  flows 
(e.g.  wakes);  flow  pattern  selection  in  flow-structure  interaction  problems;  and  vorticity 
generation  in  free  surface  flows. 
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3.3  NEW  ALGORITHMS 

Several  new  algorithms  were  developed,  including  microfluidic  meshless  techniques,  the 
boundary  cloud  method,  and  reduced  order  models  for  fast  dynamic  analysis  of 
microfluidics.  Other  new  algorithms  included  a  microflow  model  for  the  molecular 
analysis  of  microfluidic  devices  and  multiscale  methods  for  efficient  analysis  of 
microfluidic  devices. 

The  new  generation  of  spectral/hp  methods  utilizes  unstructured  meshes,  consisting  of  a 
mix  of  triangles/quadrilaterals  in  two-dimensions  or  tetrahedra/pyramids/hexahedra  in 
three-dimensions,  allowed  for  efficient  dynamic  re-discretization  "on-the-fly"  without  the 
frequent  need  for  costly  re-meshing.  An  initial  mesh  macro-skeleton  consisting  of  many 
subdomains  can  remain  unaltered  while  spectral  (p-type)  refinement  takes  place  to 
accommodate  new  features  in  the  solution.  It  is  only  when  very  large  deformations  occur 
that  new  subdomains  need  to  be  introduced,  and  given  the  unstructured  character  of  the 
discretization,  this  h— type  refinement  can  be  local  only,  preventing  global  propagation 
which  will  greatly  disrupt  load  balancing. 


3.3.1  Microfluidic  Meshless  Techniques 

Meshless  techniques  were  developed  for  the  efficient  analysis  of  microfluidics  and 
MEMS  applications  when  microfluidics  is  one  of  the  energy  domains.  The  routine 
method  of  making  a  mesh  for  a  simulation  can  be  time  consuming,  particularly  when  the 
geometry  of  the  object  to  be  modeled  is  complex. 


Figure  6.  Mesh  (a)  and  meshless  (b)  approaches 


Meshless  meshes  are  a  radically  simpler  approach  to  MEMS  modeling  and  design. 
Meshing  requirements  are  very  complemented  for  multiphysics  or  coupled  energy 
domains,  where  the  meshing  demands  are  very  acute.  The  existing  approach  of  using  a 
mesh  requires  discretization  of  the  domains  and  element  structure  (Figure  6a).  With 
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MEMS  modeling,  there  are  moving  boundary  and  strueture  problems.  An  alternate 
method  is  to  use  a  mesh  free  approaeh.  This  new  meshless  approach  requires  only 
points  (6)  or  particles,  and  there  is  no  connectivity  between  the  particles  (Figure  6b).  The 
meshless  approach  will  result  in  highly  accurate  solutions  based  upon  adaptivity  and  error 
estimation,  and  is  highly  efficient  to  implement  on  parallel  computers. 

3.3.1. 1  Finite  Point  Method 

A  cross  channel  is  one  of  the  fundamental  geometries  considered  for  the  lab-on-a-chip 
concept.  The  electroosmotic  fluidic  transport  (7)  though  a  cross-channel  was  simulated 
using  a  meshless  code  based  on  a  weighted  least  square  interpolation  is  presented  for  the 
simulation  of  electroosmotic  transport  in  capillaries.  The  geometry  of  the  cross-channel 
is  shown  in  Figure  7.  The  point  distribution  employed  for  the  cross-channel  and  the 
velocity  vector  plot  are  presented  in  Figures  8  and  9,  respectively. 


H  =  50h 


Figure  7.  The  geometry  of  the  cross-channel 
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Figure  8.  The  point  distribution  for  the  eross-ehannel 
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Figure  9.  The  resultant  veloeity  veetors  for  symmetrie  applied  potential  for  the  eross- 

ehannel 
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These  results  matched  analytical  results  published  previously  for  a  cross-shaped 
channel.  The  advantage  of  this  method  over  other  PDE  solvers  is  the  computational 
savings  due  to  a  lack  of  mesh  generation  and  integration. 

3.3.1.2  Finite  Cloud  Method 

A  finite  cloud  method  was  developed  based  upon  the  meshless  approach,  which  was 
based  on  a  fixed-kernel  technique  (8)  for  construction  of  interpolation  functions  and  a 
collocation  technique  for  discretization  of  the  Stokes  and  Navier-Stokes  equations.  The 
finite  cloud  (9)  method  uses  a  point  distribution,  and  constructs  interpolation  functions 
without  assuming  any  connectivity  between  the  points.  A  collocation  approach  is 
employed  to  obtain  a  solution  at  every  point  within  the  domain. 
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Figure  10.  Finite  cloud  injection  flow  velocity  vectors 

The  finite  cloud  method  was  applied  to  the  analysis  of  several  microfluidic  devices, 
resulting  in  radically  simpler  analysis.  For  example,  in  Figure  10  is  presented  the  results 
of  the  injection  velocity  vectors  determined  using  the  meshless  finite  cloud  method.  In 
Figure  1 1  are  presented  the  finite  cloud  method  results  of  the  horizontal  velocity  (a)  and 
the  corresponding  pressure  drop  (b)  for  this  same  injection  flow. 
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Figure  11.  Finite  cloud  results  for  horizontal  velocity  (a)  and  pressure  drop  (b)  for  an 

injection  flow 

In  Figure  12  is  presented  another  application  of  the  finite  cloud  method.  In  this  case  a 
piezoelectric  device  was  modeled  (10).  The  governing  equations  for  the  piezoelectric 
devices  were  coupled  electrical  and  mechanical  equations,  and  the  results  obtained  were 
for  the  application  of  1000  V  on  the  electrodes.  The  piezoelectric  material  was  polarized 
in  the  vertical  direction,  and  the  vertical  axial  stresses  are  displayed  in  the  figure. 
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Figure  12.  Finite  cloud  results  for  the  vertical  direction  of  a  piezoelectric  device  before 
polarization  (a)  and  after  polarization  (b)  for  an  applied  bias  of  1000  volts 


3.3. 1.3  Boundary  Cloud  Method 

A  boundary  cloud  method  is  a  combined  scattered  point/boundary  integral  approach  for 
the  boundary  only  analysis  of  exterior  problems  (11).  The  boundary  cloud  method  has 
been  combined  with  a  singular  value  decomposition  based  on  fast  algorithm  for  the  rapid 
solution  of  the  linear  system.  This  fast  boundary  cloud  method  (12)  can  be  used  for  rapid 
analysis  of  electrostatic  (13)  and  Stokes  equations.  It  was  demonstrated  the  analysis  of 
electrostatic  equations. 


3.3.2  Reduced-Order  Models  for  Fast  Dynamic  Analysis  of  Microfluidics 

A  new  reduced-order  model  for  fast  dynamic  analysis  of  microfluidics  was  developed, 
based  upon  a  weighted  Karhunen-Loeve  decomposition  approach  for  efficient  and 
accurate  dynamic  analysis  of  electro-osmotic  flows.  The  new  approach  requires  much 
fewer  number  of  modes  compared  to  the  classical  reduced-order  models,  and 
dramatically  reduced  the  central  processor  unit  (cpu)  time  required  to  design  highly 
complex  microfluidic  devices. 
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3.3.3  Microflow  Model  for  the  Molecular  Analysis  of  Microfluidic  Devices 

A  new  microflow  model  for  the  molecular  analysis  of  microfluidic  devices  was 
developed.  This  microflow  model  used  the  Direct  Simulation  Monte  Carlo  (DSMC) 
method  for  the  fundamental  physical  analysis  of  fluidic  flows  through  narrow  charmels 
(14).  Specifically,  continuum  theories  can  break  down  gas  flows  through  narrow  regions. 
Rarefaction  effects,  such  as  velocity  slip  and  temperature  jumps  that  have  been 
microscopically  and  experimentally  observed,  were  modeled  using  DSMC. 

Fluid  flows  through  microfdters  were  studied  by  DSMC.  Microfdters  designs  are  used  to 
capture  airborne  particles  such  as  bacteria,  spores  etc.  The  geometry  of  the  simulated 
fdters  is  shown  in  Figure  13.  Three  different  sized  fdters  were  analyzed,  and  the 
dimensions  of  those  filters  are  presented  in  Table  I. 

The  simulations  were  started  from  a  vacuum  initial  condition.  The  number  of  iterations 
was  determined  so  that  after  the  initial  shock  transversed  the  simulation  domain,  enough 
iterations  had  occurred  to  take  noise  free  averages.  All  simulations  were  done  on  the  NT- 
Supercluster  and  the  Origin  2000  system. 
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Figure  13.  Geometry  of  the  simulated  filter  element 
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Table  I.  Dimensions  for  Filter  Elements  (microns) 
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The  density,  pressure,  y-velocity  component  and  the  x-velocity  component  of  the  velocity 
for  a  microfilter  of  similar  geometry  are  shown  in  Figures  14-17,  respectively.  Filter 
elements  were  simulated  for  a  pressure  difference  of  0.15  Atm.  In  these  figures  the  axes 
are  in  meters  and  the  units  of  color  axes  are  noted  in  the  captions. 

In  Figures  18  and  19  are  presented  the  changes  in  pressure  and  velocity  along  the  flow 
direction  though  the  centerline  of  the  simulation  domain.  The  individual  curves  are 
shifted  along  the  x-axis  to  align  the  centers  of  the  channels. 

It  was  observed  that  for  all  cases  the  pressure  drop  was  very  close  to  linear.  This  was 
consistent  with  the  linear  increase  in  flow  velocity  along  the  channel  observed  for  filter 
elements  2  and  3.  For  filter  element  1,  the  velocity  profile,  together  with  the  linear 
density  drop,  indicated  a  degree  of  nonlinearity.  DSMC  is  an  important  technique  to 
provide  fundamental  insight  into  the  operation  of  the  microfilters. 
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Figure  14.  Density  variations  for  a  typieal  mierofilter 


Figure  15.  Pressure  variations  for  a  typieal  mierofilter  (kPa) 
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Figure  16.  Y- velocity  variations  for  a  typical  microfilter  (m/s) 
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Figure  17.  X- velocity  variations  for  a  typical  micorfilter  (m/s) 
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Figure  18.  Variation  of  pressure  along  the  mid-line 


Figure  19.  Variation  of  velocity  along  the  mid-line 
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3.3.4  Multiscale  Methods  for  Efficient  Analysis  of  Microfluidic  Devices 

Multiscale  methods  for  the  efficient  analysis  of  microfluidic  devices  have  been 
developed,  which  are  critical  for  the  efficient  and  accurate  analysis  of  microfluidic 
devices.  A  multiscale  method  is  an  approach  where  continuum  Stokes  and  Navier-Stokes 
equations  are  combined  with  molecular  approaches  such  as  Direction  Simulation  Monte 
Carlo  (DSMC)  methods  (15). 

It  was  observed  in  the  studies  of  microfilters  that  continuum  theories  broke  down  only  in 
certain  critical  regions,  and  continuum  theories  were  valid  everywhere  else.  Complex 
microfilter  designs  were  simulated  with  multiscale  approaches.  The  design  of  dual-stage 
filter  arrays  for  filtering  particles  of  a  certain  size  is  practically  impossible  by  using 
purely  DSMC  methods.  By  using  multiscale  methods,  it  was  possible  to  perform 
fundamental  studies.  These  multiscale  methods  were  at  least  100  times  faster  than  the 
purely  molecular  approaches  based  on  DSMC. 

Modeling  by  coupling  continuum  and  atomic  scale  algorithms  for  a  simple  Couette  flow 
example  is  presented.  In  Couette  flow  (Figure  20)  the  fluid  flows  between  two  walls,  one 
of  which  is  fixed,  while  the  second  wall  is  moving.  The  exact  solution  in  this  case  is 
linear,  and  differences  due  to  incompatibility  in  the  interface  conditions  cause  the 
variation  from  linearity.  The  results  of  this  DSMC  analysis  are  presented  in  Figures  21- 
23.  The  green  lines  in  these  figures  show  the  exact  solution  to  the  Couette  flow,  while 
the  DSMC  simulation  solutions  are  the  red  dots.  The  solution  at  the  interface  matches 
with  the  DSMC  solution,  meaning  that  the  interface  conditions  for  the  modeling  are 
working  properly  and  accurately. 


Wwli 

i>sivic: 

solution 

♦  . 

-  Wall 

Figure  20.  Diagram  of  Couette  flow 
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Figure  21.  DSMC  velocity  results  of  Couette  flow 


Figure  23.  DSMC  temperature  results  of  Couette  flow 
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4. 


COUPLED  CIRCUIT/DEVICE  MODELING 


The  overall  goal  of  coupled  circuit/device  modeling  was  the  development  of  a 
coupled  circuit/microfluidic  device  simulator.  Because  of  the  tight  coupling  between  the 
dynamics  of  the  fluid  flow  and  the  control  circuit,  the  response  of  the  complete  system 
can  be  accurately  determined  only  from  a  solution  of  the  circuit  equations  coupled  to  the 
microfluidics  solver  (numerical  model).  Such  a  simulation  is  unique  in  that  it  allows 
accurate  simulation  of  a  complete  microfluidic  system,  including  the  associated 
electronics,  whereby  design  tradeoffs  can  be  readily  evaluated.  Summaries  of  coupled 
circuit/device  modeling,  and  more  details  can  be  found  in  references  19-29. 


4.1  PROTOTYPE  COUPLED  CIRCUIT  AND  MICROFLUIDIC 
SIMULATORS 

Prototype  coupled  circuit  and  microfluidic  simulators  were  developed  based  on  the 
circuit  simulator  SPICES  and  the  fluidic  simulators  NsKxar  and  EOFLOW.  The 
capabilities  of  this  coupled  simulator  for  a  simple  example  illustrated  by  the  block 
diagram  in  Figure  24.  This  system  is  made  up  of  a  micropump,  a  microflow  sensor  and 
an  electronic  control  circuit.  The  electronic  circuit  adjusts  the  pump  flow  rate  so  that  a 
constant  flow  rate  is  maintained  in  the  channel.  The  flow  sensor  senses  the  flow  rate  that 
is  controlled  by  the  electronic  circuit  controlling  the  pump. 


Fluid  in 


Fluid  out 


Figure  24.  Block  diagram  of  a  microfluidic  system 


For  the  purposes  of  this  example  it  is  assumed  that  the  micropump  has  a  piezoelectric 
actuation  for  its  membrane  and  the  microflow  sensor  is  an  anemometer  type  of  flow 
sensor.  There  are  four  different  physical  domains  (electrical,  structure,  fluid,  and 
thermal)  that  must  be  considered  for  a  complete  system  level  simulation.  These  domains 
are  coupled  to  one  another  as  follows: 

•  Electromechanical  coupling  for  piezoelectric  actuation  of  the  pump  membrane 

•  Fluid-structure  coupling  due  to  volume  displacement  of  the  pump  membrane 

•  Fluid-thermal  coupling  because  of  the  thermoresistor  cooling  the  fluid  for  an 
anemometer  type  of  microflow  sensor 

•  Electrothermal  coupling  for  thermoresistor  heating  due  to  current  flow  in  the 
microflow  sensor 
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The  overall  system  is  simulated  using  a  eombination  of  eoupled  solvers,  compact 
models  and  lumped  elements.  Simple  elements  are  implemented  as  lumped  elements 
or  compact  models.  These  elements  are  the  electromechanical  transducer  (piezoelectric 
actuator)  and  thermoresistors. 


4.2  PROTOTYPE  FLOW  SENSORS 

Flow  sensors  are  much  more  complicated  but  often  the  fluid  flow  in  sensors  is  relatively 
simple.  For  example,  if  the  fluid  flow  has  a  constant  parabolic  profile  for  the  velocity, 
then  it  is  possible  to  use  compact  models  for  the  flow  sensors  as  well.  It  is  important  to 
note  that  these  compact  models  are  parameterized  and  can  be  highly  nonlinear.  These 
models  are  obtained  by  insight  gained  from  detailed  physical  level  simulations.  Usually 
pumps  operate  in  a  nonlinear  mode  of  fluid  flow  with  a  strong  fluid- structure  interaction. 
Therefore,  a  detailed  physical  level  simulation  of  the  pump  is  required. 

In  this  example,  the  electromechanical  actuator,  thermoresistors  and  flow  sensors  are 
described  as  lumped  elements  and/or  compact  models.  The  pump  is  modeled  at  the 
detailed  physical  level.  All  lumped  elements  and  models  are  implemented  in  SPICE3. 
The  pump  is  implemented  as  a  direct  SPICES-NsKxar  interconnection  (Figure  25). 
SPICE3  transfers  the  time  tspice  and  pressure  P  for  membrane  actuation  and  to  receive  the 
flowrate  Q  and  the  time  tcaii  for  the  next  call  from  NsKxar. 


PUMP  model 


Figure  25.  The  SPICE3-  interaction  for  the  pump  microsystem 
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The  model  for  eleetromechanieal  eoupling  with  a  piezoeleetrie  actuation  of  the 
membrane  is  shown  in  Figure  26.  This  model  forms  the  interface  between  the 
electrical  and  mechanical  networks.  The  electrical  characteristics  of  the  piezoelectric 
actuator  are  described  by  the  capacitor  C.  The  input  voltage  v  translates  into  an  output 
pressure  P  by  virtue  of  the  piezoelectric  effect  with  coefficient  k  to  activate  the  pump 
membrane.  The  mechanical  characteristics  of  the  piezoelectric  actuator  are  coupled  with 
the  mechanical  characteristics  of  the  substrate. 


Figure  26.  The  lumped  model  for  piezoelectric  actuation 


4.3  ANEMOMETER  FLOW  SENSOR 

An  anemometer  type  microflow  sensor  Macromodel  has  been  developed  using  neural 
networks.  The  numerical  solutions  of  Partial  Differential  Equations  (PDE’s)  describing 
the  microflow  sensor  were  used  as  training  data  for  the  neural  networks.  This  model  is 
simple  and  accurate  and  has  been  verified  with  numerical  simulations  for  a  wide  range  of 
sensor  parameters.  The  model  incorporates  both  the  steady  state  and  the  dynamic 
responses  of  the  flow  sensors.  The  model  has  been  implemented  in  SPICES  and  can  be 
used  for  simulating  the  flow. 

In  Figure  27  are  presented  PDE  temperature  solutions  for  a  microflow  sensor  with  zero 
fluid  flow  in  the  channel  (a),  and  with  flow  (b).  The  height  of  the  channel  was  20 
microns.  In  Figure  28  are  presented  the  PDE’s  modeling  results  of  the  change  in 
temperature  verses  distance  for  a  50-micron  channel  height.  Superimposed  on  this  figure 
is  experimentally  measured  temperature  differences,  which  validates  the  PDE  model. 
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Figure  27.  PDF  resultant  temperatures  (Kelvin)  for  a  microflow  sensor  with  zero  flow 

(left)  and  with  flow  (right) 


d  (|u.m) 


Figure  28.  The  prediction  accuracy  of  the  neural  network  based  model  for  a  channel 

height  of  50  pm 
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4.4  COMPLETE  SYSTEM  FOR  SIMULATION 


The  structure  of  a  complete  anemometer  type  flow  sensor  is  shown  in  Figure.29,  a 
macromodel  has  been  developed  (16,17).  This  macromodel  is  based  on  neural  networks 
trained  using  data  (Figure  30)  from  detailed  physical  simulations. 


Figure  29.  The  structure  of  an  anemometer  type  microflow  sensor 
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Steady  state  nominal  model  dynamic  extension 

Figure  30.  The  dynamic  macromodel  for  the  flow  sensor 


The  steady-state  solution  Tsso  corresponds  to  a  nominal  power  for  the  heat  source  %.  The 
neural  network  output  Tsso  is  a  multivariate  function  of  the  flow  velocity  U  and  the 
vector  of  geometrical  and  physical  parameters  0.  Tss  is  a  linear  function  of  the  heat 
source  X  and  Tsso- 

The  inputs  to  the  neural  network  are  the  flow  velocity  U  and  the  vector  of  geometrical 
and  physical  parameters  ©.  The  results  from  this  model  are  in  good  agreement  with  the 
simulated  data  for  a  large  range  of  parameters.  The  dynamic  macromodel  is  incorporated 
in  SPICES  by  coupling  with  a  sensor  circuit  and  a  model  for  thermoresistors  for  the 
heater  and  sensors. 

A  simplified  schematic  of  the  complete  system  for  simulation  is  shown  in  Figure  31.  In 
this  system,  the  flow  rate  Q  determines  the  flow  sensor  velocity  U  for  a  given  set  of 
geometry  parameters  (h,  d,  wsens).  Based  on  the  fluid  flow  rate  the  thermoresistor 
temperatures,  Tl,  T2  and  T3,  change  which  in  turn  alters  the  resistance  values  Rl(Tl), 
R2(T2)  and  R3(T3).  The  resistance  Rl(Tl)  and  R3(T3)  are  included  in  a  Wheatstone- 
bridge  arrangement  with  two  fixed  resistors  R4  and  R5.  The  voltage  difference,  Vr3-Vri, 
is  directly  proportional  to  the  temperature  difference  T3-T1.  This  voltage  difference  is 
linearly  transformed  to  the  output  voltage  Vout  by  an  operational  amplifier  with  a 
controlled  gain.  This  output  voltage  determines  the  pressure  P,  which  activates  the  pump 
membrane  and  changes  the  flow  rate  Q.  The  thermoresistor  of  the  heater  (R2)  is 
activated  by  the  control  electronics  that  maintains  a  constant  heater  temperature. 
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Figure  31.  A  description  of  the  complete  system  for  simulation 

The  circuit  simulator  SPICES  is  chosen  as  the  controlling  solver  since  it  has  automatic 
time  step  control,  and  supports  lumped-element  equivalent  circuits.  Furthermore,  models 
for  different  abstraction  levels  can  be  easily  implemented  in  SPICE3.  is  embedded  as  a 
subroutine  in  SPICES.  The  interaction  with  SPICES  is  by  means  of  the  model  code  and 
the  simulation  engine.  Synchronization  timepoints  are  determined  and  used  by  the 
SPICES  transient  analysis  engine.  The  pump  is  modeled  as  a  SPICES  element  with  being 
the  underlying  simulation  engine.  Lumped  element  descriptions  and/or  compact  models 
describe  the  other  elements  in  the  circuit. 

The  simulation  results  from  the  coupled  simulator  are  presented  in  Figure  32.  In  this 
simulation,  one  can  determine  the  pressure  on  the  pump  membrane,  the  flow  velocity, 
and  the  output  control  voltage  as  a  function  of  time  for  various  component  parameters.  In 
this  simulation,  the  structure  of  the  pump  and  the  structure’s  influence  on  overall  system 
performance  can  be  easily  determined.  The  coupled  simulator  provides  valuable 
information  to  the  system  or  device  developer. 

The  flow  pump  rate  Q  determines  the  flow  sensor  velocity  U.  This  yields  the 
temperatures  for  the  sensor  thermoresistors.  The  difference  between  the  resistance  values 
Rl(Tl)  and  R3(T3)  is  transformed  into  the  voltage  Vout  by  the  control  electronics.  Vout 
is  used  to  control  the  pressure  P  for  the  pump  membrane  that  in  turn  determines  the  flow 
rate  Q. 
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Figure  32.  The  pressure,  velocity  and  output  voltage  for  the  simulation  example  as  a 

function  of  time 


SPICE3  and  EOFLOW  simulator  codes  were  modified  to  enable  coupled  simulation  with 
any  number  of  interconnected  cross  channels  (Figure  33).  The  results  for 
SPICE3/EOFLOW  coupling  are  briefly  described.  A  system  described  by  a  mesh  of 
10x10  interconnected  cross  channel  pipes  was  simulated  to  demonstrate  the  advantages 
of  the  coupled  simulation  environment.  Voltage  sources  with  different  amplitude  and 
frequency  were  connected  to  every  outlet  as  shown  in  Fig.  y.  The  response  was  checked 
near  the  center  of  the  mesh.  A  transient  analysis  for  300  ps  took  an  hour  of  simulation 
time  on  a  Sun  Solaris  with  UltraSparc  440  MHz  processor.  The  goal  of  this  simulation 
was  to  check  the  voltage  calculation  at  the  interconnection  point  of  the  two  channels. 
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The  resulting  output  signal  eontains  all  harmonies  eorresponding  to  input 
frequeneies  with  amplitudes  as  expeeted.  The  goal  of  this  simulation  was  to  eheek 
the  voltage  ealeulation  at  the  intereonneetion  point  of  the  two  ehannels.  The  resulting 
output  signal  eontains  all  harmonies  eorresponding  to  input  frequeneies  with  amplitudes 
as  expeeted.  Figure  34  shows  one  period  of  output  waveform  and  Figure  35  eontains  the 
speetrum  of  output  waveform. 


1000*sin(2*pi*100000)V 

900*sm(2*pi*90000)V 

800*sin(2*pi*80000)V 

700*sin(2*pi*70000)V 

600*sin(2*pi*60000)V 

50Q*sin(2*pi=^50000)V 

40a*sin(2*pi*40000)V 

30a*sin(2*pi*30000)V 

200^sin(2*pi*20000)V 

10a*sin(2*pi*i0000)V 


Figure  33.  10x10  eross  ehannel  mesh  simulated  with  applied  voltages  at  all  outlets 


29 


2 


4  6 

Frequency  (Hz) 


Figure  35.  Spectrum  of  the  output  signal 


SPICE3/E0FL0W  simulator  codes  were  modified  to  implement  a  new  crosschannel 
model  in  SPICE.  Other  code  modifications  were  made  to  enable  coupled  simulation  with 
two  interconnected  crosschannels.  A  microfluidic  system  containing  two  cross  channel 
pipes  (Figure  36)  was  simulated.  The  goal  was  to  check  the  voltage  calculation  at  the 
interconnection  point  of  the  two  channels.  These  results  match  the  previous  results 
obtained  for  the  simulation  of  one  cross  channel  pipe. 
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Figure  36.  Example  of  two  eross  ehannel  pipes  simulated  during  one  of  the  time  steps  in 

the  transient  analysis 
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4.5  ELECTROOSMOTIC  FLOW  SIMULATOR 


The  code  for  the  electroosmotic  flow  simulator  was  modified  to  correct  memory  leakage 
problems.  Other  code  modifications  were  made  to  provide  reliable  and  efficient 
circuit/microfluidic  device  coupling.  The  electroosmotic  flow  simulator  was  adapted  for 
repeated  calls  with  restart  capabilities.  The  circuit/microdevice  is  based  on  a  generic 
microfluidic  device  model  that  allows  use  of  different  solvers  in  a  similar  manner. 

Applications  to  microfluidic  system  design  based  on  the  electroosmotic  effect  were 
demonstrated.  A  microfluidic  system  consisting  of  a  cross  channel  pipe  driven  by  control 
electronics  (Figure  37)  was  simulated.  The  control  electronics  implements  an  integral 
type  controller  and  produced  an  output  that  is  an  integral  of  the  flow  rate  through  the 
vertical  channels.  The  goal  of  this  system  was  to  obtain  zero  flow  rate  through  the 
vertical  channels  by  automatically  determining  the  voltages  applied  to  the  vertical 
channels.  Simulation  results  (Figure  38  and  39)  show  the  dynamic  behavior  of  the  whole 
system.  This  simulation  provides  the  basis  for  future  simulation  of  separations  and 
dosing  systems  for  biomedical  and  chemical  application  such  as  a  micro  total  analysis 
system. 


Figure  37.  Example  of  a  controlled  electroosmotic  system 
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Figure  38.  Dynamic  behavior  of  the  control  potential  and  the  flow  rate  for  the 

electroosmotic  system 


Figure  39.  Dynamic  behavior  of  streamlines  for  the  electroosmotic  system 
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5.  EXPERIMENTAL  EFEORT 

Summaries  of  the  experimental  work  are  presented  here,  and  more  details  can  be  found  in 
references  30-33. 

5.1  FABRICATION 

New  materials  have  been  used  for  the  first  time  in  MEMS  fabrications,  and  novel  devices 
have  been  fabricated  for  the  first  time. 

5.1.1  Materials 

Liquid  crystals  and  biodegradable  polymers  have  been  used  as  microfluidic  materials. 

5.1. 1.1  Liquid  Crystal  Polymer  as  a  Microfluidics  Material 

Recently,  new  techniques  using  low-cost  organic  polymers  have  been  studied.  Fluid 
circuits  realized  using  photoresist,  parylene,  and  PDMS  (18)  have  been  demonstrated. 
However,  a  common  difficulty  encountered  in  preparing  PDMS  microfluidic  devices  is 
the  relative  weak  bonding  between  PDMS  and  other  materials  -  lateral  leakage  is  a 
major  concern  for  such  fluid  devices. 

The  liquid  crystal  polymer  (LCP)  is  a  thermoplastic  polymer  material  with  unique 
structural  and  physical  properties.  It  contains  rigid  and  flexible  monomers  that  link  to 
each  other.  When  flowing  in  liquid  crystal  state,  the  rigid  segments  of  the  molecules 
align  next  to  one  another  in  the  direction  of  shear  flow.  Once  this  orientation  is  formed, 
their  direction  and  structure  persist,  even  when  LCP  is  cooled  below  the  melting 
temperature.  This  is  different  from  most  thermoplastic  polymers  (for  example  Kapton), 
whose  molecules  are  randomly  oriented  in  the  solid-state. 

An  alternative  process  for  fabrication  a  fluid  channel  using  LCP/glass  assembly  was 
developed  (19).  Compared  with  other  methods  of  forming  glass  fluid  channels,  LCP 
has  several  unique  advantages.  First,  LCP  fdms  are  virtually  impermeable  by  liquid 
and  gas.  Such  a  property  may  be  important  in  situations  where  liquid  and  gas  diffusion 
is  undesired.  Secondly,  whereas  PDMS  can  be  bond  to  glass  by  natural  adhesion  force, 
LCP  can  be  attached  to  glass  with  much  stronger  bonding  developed  through  the 
lamination  process.  Such  bonding  strength  is  important  for  retaining  the  integrity  of  the 
fluid  chip.  Commercial  LCP  fdms  are  also  with  a  lower  cost  compared  to  PDMS. 
Thirdly,  LCP  to  glass  bonding  requires  much  lower  temperature  and  shorter  processing 
time  compared  with  glass-glass  bonding. 

5.1. 1.2  Pyrex  Glass  Wafers  as  a  Microfluidics  Material 

Microchannels  on  a  pyrex  glass  wafer  were  fabricated  for  the  first  time  (20).  This  glass 
wafer  was  patterned  by  using  photolithography  process  with  gold  mask  and  etched  in 
HF  to  form  a  70-//m  deep  trench.  A  LCP  film  was  then  thermally  laminated  on  the 
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glass  wafer  with  the  glass  heated  to  260°C  on  a  hotplate.  The  lamination  of  TCP  is 
conducted  by  applying  moderate  pressure  using  a  roller.  Via  holes  through  LCP 
fdm  to  glass  channel  reservoirs  were  produced  using  mechanical  punching  or  laser 
drilling.  Plastic  inlet  and  outlet  pipes  were  bonded  on  LCP  fdm  with  epoxy  as  the 
reservoirs. 

The  functionality  of  such  micro  fluid  channels  in  electrophoresis  tests  was 
demonstrated  (21).  The  LCP-glass  assembly  was  first  immersed  in  TBE  buffer  solution 
and  degassed  20  minutes  in  a  vacuum  chamber  to  fill  the  entire  channel  and  remove 
trapped  air  bubbles  (22).  Fluorescent  beads  were  then  injected  into  the  fluid  channel 
through  one  of  the  inlets.  Platinum  electrodes  were  inserted  at  the  inlet  and  outlet 
reservoirs.  The  hydraulic  potential  was  completely  balanced  before  electrophoresis  was 
processed.  When  50  V/cm  was  applied  across  the  channel,  the  speeds  of  flowing 
fluorescent  beads  varied  from  36-60  pm/s  at  different  locations  in  the  channel  system. 

As  an  alternative  microchannels  in  LCP  film  were  created.  An  Aluminum  layer  was  first 
deposited  on  LCP  film.  This  A1  layer  was  patterned  by  photolithography  to  form  a  mask 
of  micro  channel.  The  LCP  film  was  then  etched  to  desired  shape  and  depth  by  RIE 
using  the  A1  mask.  After  etch  completion,  the  A1  mask  was  stripped  away  to  release  LCP 
film.  This  LCP  film  with  the  fluid  channel  on  it  was  thermally  bonded  to  a  glass 
substrate.  An  SEM  micrograph  (Figure  40)  clearly  demonstrates  the  cross-section  view 
of  the  channel.  This  potentially  allows  LCP  channels  to  be  bonded  to  other  flexible 
films,  such  as  another  LCP  film. 
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Figure  40.  (a)  SEM  micrograph  of  a  sealed  LCP  microchannel;  (b)  images  of 
fluorescence  beads  moving  inside  a  micro  channel 


5.1. 1.3  Biodegradable  Polymer  as  a  Microfluidics  Material 

We  report  the  development  of  micromachining  techniques  for  a  biodegradable  polymer 
for  the  first  time  (23,24).  By  virtue  of  their  ability  to  naturally  degrade  in  tissues, 
biodegradable  polymers  hold  immense  promise  as  new  materials  for  implantable 
biomedical  microdevices.  This  work  focuses  on  establishment  of  microfabrication 
processes  for  biodegradable  microstructures  and  microdevices.  Three  unique  fabrication 
processes  have  been  established:  (1)  micro-molding  process  to  form  3D  microstructures 
in  polycaprolactone  (PCL)  via  a  silicon  micromachined  mold;  (2)  a  method  of 
transferring  metal  patterns  to  surfaces  of  PCL  substrates;  (3)  techniques  for  sealing  both 
dry  and  liquid-filled  PCL  micro-cavities  with  a  metal  thin  film  (e.g.  gold).  Chemical 
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compatibility  of  PCL  with  common  micromachining  chemicals  has  been 
investigated  (Figure  41). 


[c] 


Figure  41.  (a)  SEM  mierograph  of  a  mieromachined  channel  in  BDP;  (b)  fabrieation 

proeess  for  a  sealed  cavity  in  BDP 
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5.1.2  Microfluidic  Components 


Microfluidic  components  including  channels,  pumps  and  valves  were  fabricated.  Flow 
visualization  methods  were  developed  to  characterize  these  components.  Also,  a  micro 
fabricated,  externally  powered  magnetic  diaphragm  pump  for  telemetric  operation  of 
microfluidic  systems  was  fabricated,  as  was  a  micro  magnetic  mixer  for  rapid,  on- 
demand  mixing  of  fluidics. 


5.1.2.1  Magnetic  Micro  Pump 

A  tetherless  micro  magnetic  membrane  pump  was  fabricated  capable  of  remote  operation 
using  an  external  magnetic  field  (due  to  the  long-range  forces  provided)  (25-28).  The 
basic  structure  of  the  fully  assembled  pump  is  illustrated  in  Figure  42,  which  shows 
Layer  I  placed  on  top  of  Layer  If.  A  fluid  flow  rate  of  1.2-pl/min  was  measured  at  a 
pumping  frequency  of  2.9-Hz  in  the  presence  of  a  2.85xlO^-A/m  oscillating  external 
magnetic  field. 


Figure  42.  Schematic  diagram  of  a  magnetic  actuated  micro  pump 


Layer  I,  which  consisted  of  the  diffuser  elements  (which  control  one-way  fluid  flow), 
pump  chamber,  micro  channels,  and  inlet/outlet  wells,  was  fabricated  using  polymer 
micromachining  techniques.  The  material  used  is  a  flexible,  transparent  silicone 
elastomer,  polydimethyl  siloxane  (PDMS).  Due  to  strong  PDMS-PDMS  bonding.  Layer 

1  can  thus  be  self-bonded  to  Layer  11  (which  has  a  thin  PDMS  membrane)  by  pressing  the 

2  layers  together. 

Layer  11,  the  magnetic  membrane  actuator,  is  made  of  ferromagnetic  materials  embedded 
within  a  thin  PDMS  membrane.  In  the  presence  of  the  external  magnetic  field,  the 
generated  magnetic  torque  causes  the  ferromagnetic  pieces  to  deflect,  thereby  displacing 
the  membrane  as  shown  in  Figure  42.  Maximum  membrane  displacement  measured 
(using  our  present  strongest  magnet)  is  83.9  pm  under  a  2.85xlO^-A/m  external  magnetic 
field.  Membrane  displacement  thus  pushes  fluid  out  of  the  pump  chamber,  and  an 
oscillating  magnetic  field  can  be  used  to  affect  pumping  action. 
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This  pump  is  suitable  for  implantable  biomedical  microfluidic  systems  as  (1)  PDMS 
is  physically  and  chemically  stable,  thus  it  is  safe  and  more  compatible  than  silicon 
for  an  implanted  device,  (2)  material  damage  is  only  of  concern  for  ~  5x  higher 
magnitudes  of  applied  magnetic  fields  than  the  ones  we  used  for  our  device,  and  (3)  the 
capability  for  remote  operation  (without  the  need  for  tether  wires  for  power  transmission 
and  pump  actuation),  when  combined  with  a  wireless  communication  module,  realizes  a 
truly  remotely  activated  drug  delivery  chip. 

5.1.2.2  Magnetic  Micro  Mixer 

A  micromachined  magnetic  bar  stirrer  has  been  developed,  and  the  basic  schematic  of  the 
bar  stirrer  is  shown  in  Figure  43.  The  micromagnetic  bar  is  anchored  by  a  hub.  It  is 
magnetized  by  an  external  magnetic  field,  and  is  capable  of  rotating  under  a  remotely 
applied  magnetic  field.  It  was  inspired  by  a  macroscopic  bar  mixer. 

Most  microfluidic  mixers  are  coupled  with  a  flow  channel,  but  not  with  a  fluid  reservoir 
or  a  reaction  chamber.  As  a  result,  common  biological  procedures  such  as  the  mixing  of 
two  components  are  impossible.  The  magnetic  bar  mixer  is  motivated  by  the  needs  of 
biologist  and  chemists  to  have  controlled  mixing. 

Currently  a  second-generation  device  is  being  developed.  Although  bar  rotation  has  been 
demonstrated,  smooth  rotation  has  not  been  accomplished.  The  roughness  of  the  mask 
contributes  to  this  shortcoming.  In  the  next  generation  device,  the  process  and  design 
will  be  optimized  based  on  existing  knowledge. 


micro  magnetic  stirrer  hub/anchor 


Figure  43.  Schematic  diagram  of  a  micromachined  micro  mixer 


In  Figure  44  is  presented  a  numerical  simulation  of  the  fluid  flow  forced  by  the  moving 
magnetic  bar  in  an  effort  to  optimize  the  mixing  performance.  The  simulation  plot 
(Figure  45)  illustrates  tow  fluids  being  mixed  with  a  rotating  magnetic  bar  stirrer. 
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Figure  44.  Illustration  of  two  fluids  being  mixed  with  a  rotating  magnetie  bar  stirrer 


Figure  45.  Sequential  images  of  two  fluids,  eolored  blue  and  yellow,  being  mixed  in  a 
mieroehannel  with  a  miero  magnetie  bar  mixer 
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5.2  ELECTROKINETIC  MICROCHANNEL  GEOMETRIES 

Electrokinetic  microchannel  geometries  have  been  optimized,  and  novel  flow 
diagnostics  have  been  developed,  and  are  discussed  in  this  chapter. 


5.2.1  Design  and  Optimization 

A  comprehensive  effort  to  design  and  optimize  electrokinetic  turn  geometries  for 
minimally  dispersive  capillary  electrophoretic  separations  was  undertaken  and 
completed.  This  analysis  included  the  development  of  an  analytical  model  and  the 
implementation  of  novel  diagnostic  techniques  to  confirm  the  predictions  from  the 
analytical  model.  Two  techniques,  bleached  fluorescence  and  caged  fluorescence,  were 
developed  and  applied  as  part  of  this  Composite  CAD  project. 

In  Figure  46  is  presented  a  caged  fluorescence  image  of  a  flow  field  with  both  pressure 
driven  and  electroosmotic  flow  components.  The  band  on  the  left  is  optically  injected 
into  the  flow  field  using  a  355  nm  laser  probe.  The  flow  of  the  fluid  is  from  left  to  right 
along  the  channel.  On  the  right  of  this  figure  is  the  same  band  at  a  later  time.  The  effect 
of  the  drag  from  the  walls  of  the  channel  is  evident  from  the  shape  of  this  band. 


Figure  46.  Caged  fluorescence  image  of  a  flow  field 

The  bleached  fluorescence  technique  (29)  was  a  novel  technique  with  advantages  that  are 
very  complimentary  to  those  of  caged  fluorescence.  Experimental  results  were  compared 
with  simulations  to  demonstrate  that  simulations  could  be  used  as  effective  design  tools 
(Figure  47). 
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t  =  0  ms 


Figure  47.  Imaging  of  the  pressure  driven  flow  in  a  microchannel  using  the  improved 
caged  fluorescence  imaging  system.  The  intensity  of  the  plots  shown  in  (b) 
correspond  to  the  centerline  of  the  first,  second  and  forth  images  shown  in  (a). 


In  an  initial  study  of  the  dispersion  problem,  an  analytical  model  describing  the 
dispersion  caused  by  constant  radius  turns  was  developed.  This  model  spanned  the 
different  regimes  that  exist  for  analyte  dispersion  in  electrokinetic  channels  in  the  context 
of  the  advective-diffusion  equation:  axial  diffusion  limit,  Taylor-Aris  limit  and  the  pure 
advection  limit.  Figure  48  shows  the  turn  variance  ratio  as  a  function  of  the  turn 
transport  ratio  for  an  180°  constant  radius  turn.  The  introduction  of  the  non-dimensional 
term,  ttum  allows  for  a  quantitative  determination  of  the  dispersion  regimes  associated 
with  the  constant  radius  turn. 
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Figure  48.  Turn  variance  ratio  verses  transport  ratio  for  a  180'’  constant  radius  turn.  The 
caged  fluorescence  images  show  the  experimentally  observed  dispersion  regimes  for 

the  different  turn  ratios 


A  novel  technique,  bleached  fluorescence  visualization,  for  imaging  microflows  was 
developed  at  the  Stanford  Microfluidics  Laboratory.  This  technique  was  demonstrated  as 
a  viable  tool  for  measuring  molecular  diffusivities  and  electrokinetic  mobilities  in  both 
glass  and  acrylic  channels  (Aclara  BioSciences).  Figure  49  shows  the  evolution  of  a 
photo-bleached  line  in  a  steady  electroosmotic  flow  field.  A  time-sequence  of  column 
averaged  axial  intensity  profiles  with  corresponding  Gaussian  fits  is  shown  in  Figure  50. 
This  method  was  demonstrated  to  yield  a  simultaneous  measurement  of  the 
electrophoretic  mobility  and  diffusion  coefficient  of  the  fluorescein  dye. 
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t  =  0  s  t  =  0.53  s 


Figure  49.  Evolution  of  a  photobleached  line  in  steady  eleetroosmotic  flow  field 


Figure  50.  Time-sequence  of  column-averaged  axial  intensity  profiles  for  the  pure 
eleetroosmotic  flow  experiment  shown  in  Figure  48 
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Both  bleached  and  caged  fluorescence  techniques  were  applied  to  validate  optimal 
electrokinetic  turns  for  both  90°  and  180°  turns.  Figure  51  shows  the  reduction  in 
dispersion  achieved  through  the  optimized  turn  geometries  compared  to  a  simple  180° 
turn.  A  significant  reduction  in  the  variance  has  been  shown  for  this  particular  geometry. 
The  image  with  the  optimized  geometry  displays  a  simultaneous  separation  of  different 
specie  present  in  the  uncaged  sample  plug  along  with  reduced  sample  dispersion. 


Figure  51.  Caged  fluorescence  (i.e.,  experimentally  obtained  visualizations)  results  for 
dispersion  of  a  sample  plug  around  the  two  cases  of  a  simple,  uncorrected  and  a 

corrected  180°  turn 


5.2.2  Design  Rules  for  Optimization 

Design  rules  for  the  optimization  of  90°  and  180°  constant  radius  turns  have  been 
proposed  based  on  simulations  and  predictions  from  an  analytical  model  developed.  This 
model  is  valid  across  all  the  dispersive  regimes  allowed  by  the  advective-diffusion 
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equation.  A  novel  experimental  technique,  bleached  fluorescence  visualization,  was 
implemented  to  verify  these  predictions  (Figure  52).  We  also  implemented  caged- 
fluorescence  visualization  in  these  experiments  and  compared  the  two  techniques.  An 
optimized  turn  geometry  shown  to  reduce  the  dispersion  variance  to  less  than  1%  of  a 
simple  turn  was  proposed  (Figure  53). 


Figure  52.  Simulations  of  microfluid  flow  in  180  degrees  microchannel  bends 


Figure  53.  Caged  fluorescence  measurements  of  flows  in  180  degree  microchannel  bends 

with  different  geometries 


46 


6. 


CONCLUSIONS 


This  project  led  to  several  breakthrough  advances  in  computational  prototyping  tools  for 
microfluidics.  Sophisticated  and  advanced  computational  design  tools  were  developed 
for  microfluidics.  These  tools  are  several  orders  of  magnitude  faster  than  existing 
simulation  tools.  In  addition,  the  radically  new  algorithms  developed  as  part  of  this 
project  have  laid  foundations  for  next  generation  design  tools  in  microfluidics.  All  the 
design  tools  have  been  rigorously  validated  with  experimental  data.  Using  the  design 
tools,  several  complex  fluidic  devices  have  been  designed  and  several  fabrication  cycles 
have  been  eliminated  through  the  use  of  these  design  tools.  The  results  of  this  effort  are 
being  transitioned  to  commercial  MEMS  CAD  companies,  and  some  of  this  work  is 
already  embedded  into  their  tools. 


47 


7.  REFERENCES 

1.  N.  R.  Alum,  “CAD  enabled  MEMS  innovation”,  Journal  of  Materials 
Processing  <&  Manufacturing  Science,  Vol.  8,  No.  4,  pp.  325-327,  April  2001. 

2.  D.  Liu,  M.  Maxey  and  G.E.  Kamiadakis,  A  fast  algorithm  for  partieulate 
mieroflows  in  eomplex  geometries  Proe.  ASME  Winter  Meeting,  NY,  November 
2001. 

3.  I.  Lomtev,  R.M.  Kirby,  and  G.E.  Kamiadakis,  "A  discontinuous  Galerkin  ALE 
method  for  compressible  viscous  flows  in  moving  domains",  J.  Comp.  Phys., 
vol.  155,  pp.  128-159,1999. 

4.  C.  Evangelinos,  D.  Lucor,  C.-H.  Su  and  G.E.  Kamiadakis,  Flow-induced 
vibrations  of  non-linear  cables.  Part  1:  Models  and  algorithms,  Int.  J.  Num. 
Meth.  Engineering,  to  appear,  2001. 

5.  D.  Lucor,  C.  Evangelinos,  L.  Imas  and  G.E.  Kamiadakis  Flow-induced 
vibrations  of  non-linear  cables.  Part  11:  Simulations  Int.  J.  Num.  Meth. 
Engineering,  to  appear,  2001. 

6.  N.  R.  Alum,  “A  point  collocation  method  based  on  reproducing  kernel 
approximations”.  International  Journal  for  Numerical  Methods  in  Engineering, 
Vol.  47,  No.  6,pp.  1083-1121,2000. 

7.  M.  Mitchell,  R.  Qiao,  and  N.  R.  Alum,  “Meshless  analysis  of  steady-state 
electroosmotic  transport”.  Journal  of  Microelectromechanical  Systems,  Vol.  9, 
No.  4,pp.  435-449,2000. 

8.  N.  R.  Alum,  “A  reproducing  kernel  particle  method  for  meshless  analysis  of 
microelectromechanical  systems”.  Computational  Mechanics,  Vol.  23,  No.  4, 
pp.  324-338,  1999. 

9.  N.  R.  Alum  and  G.  Li,  “Finite  Cloud  Method:  A  tme  meshless  technique  based 
on  a  fixed  reproducing  kernel  approximation”,  International  Journal  for 
Numerical  Methods  in  Engineering,  Vol.  50,  No.  10,  pp.  2373-2410,2001. 

10.  R.  Ohs  and  N.  R.  Alum,  “Meshless  analysis  of  piezoelectric  devices”. 
Computational  Mechanics,  Vol.  27,  No.  l,pp.  23-36,2001. 

1 1 .  G.  Li  and  N.  R.  Alum,  “Boundary  Cloud  Method:  A  combined  scattered 
point/boundary-  integral  approach  for  boundary-only  analysis”,  Accepted  for 
Publication  in  Computer  Methods  in  Applied  Mechanics  and  Engineering,  2001 . 

12.  Vaishali  Shrivastava  and  N.  R.  Alum,  “A  fast  boundary  cloud  method  for 
exterior  2-D  electrostatic  analysis”.  Submitted  for  Publication,  2001. 

13.  Gang  Li  and  N.  R.  Alum,  “Linear,  nonlinear  and  mixed-regime  analysis  of 
electrostatic  MEMS”,  Sensors  and  Actuators  A,  Vol.  91,  No.  3,  pp.  278-291, 
2001. 


48 


14.  O.  Aktas,  N.  R.  Alum,  U.  Ravaioli,  “Application  of  a  parallel  DSMC 
technique  to  prediet  flow  eharacteristics  in  microfluidie  filters”  Journal  of 
Microelectromechanical  Systems,  Vo\.  10,  No.  4,  pp.  538-549,2001. 

15.  O.  Aktas  and  N.  R.  Alum,  “A  combined  continuum/DSMC  approaeh  for 
multiscale  analysis  of  microfluidic  filters”,  Submitted  for  Publication,  2001. 

16.  R.M.  Kirby,  G.E.  Kamiadakis,  O.  Mikulehenko,  K.  Mayaram  An  Integrated 
Simulator  for  Coupled  Domain  Problems  in  MEMS,  J.  MEMS,  vol.  10  (3),  p. 
379,  2001. 

17.  R.M.  Kirby,  G.E.  Kamiadakis,  O.  Mikulehenko,  K.  Mayaram  Integrated 
Simulation  for  MEMS:  Coupling  Flow-Stmcture-Thermal-Electrical  Domains, 
Chapter  5,  The  MEMS  Handbook,  edited  by  M.  Gad-el-Hak,  CRC,  2001. 

18.  Deniz  Armani,  Chang  Liu,  "Re-configurable  Fluid  Cireuits  By  PDMS  Elastomer 
Micromachining",  12th  International  Conference  on  MEMS,  MEMS  99,  pp.222- 
227,  Orlando,  FL,  1998. 

19.  X.  Wang,  L.H.  Lu,  C.  Liu,  "Micromachining  Techniques  for  Liquid  Crystal 
Polymer",  MEMS  01,  Switzerland,  January  2001. 

20.  L.H.  Lu,  X.  Wang,  C.  Liu,  "An  efficient  low-cost  micro  fabrication  method  of 
capillary  electrophoresis  chip  using  liquid  crystal  polymer  (LCP)",  LI 006,  14th 
International  Symposium  on  Microscale  Separations  and  Analysis,  January  13-18, 
2001,  Boston,  MA. 

21. J.I.  Molho,  AE  Herr,  BP  Mosier,  JG  Santiago,  TW  Kenny,  et  al.,  “X  low 
dispersion  turn  for  miniaturized  electrophoresis’’,  Solid  State  Sensor  and 
Actuator  Workshop,  Hilton  Head  Island,  SC.  2000. 

22.  Bo  Choi,  Min  Ma,  Christopher  White,  Chang  Liu,  "Electrolytic  and  thermal 
bubble  generation  using  AC  inductive  powering,"  International  Conferenee  on 
Solid-State  Sensors  and  Actuators  (Transducer'99),  Sendai,  Japan,  7-10  June 
1999. 

23.  D.  Armani,  C.  Liu,  "Microfabrication  Technology  for  Polycaprolactone,  A 
Biodegradable  Polymer",  J.  of  Mieromechanics  and  Microengineering,  10,  pp. 
80-84,  2000. 

24.  D.  Armani,  C.  Liu,  "Microfabrication  Technology  for  Polycaprolactone,  A 
Biodegradable  Polymer"  Int.  Conf.  on  MEMS,  MEMS  2000,  Japan. 

25.  M.  Khoo,  C.  Liu,  "Micro  magnetic  silicone  elastomer  membrane  actuator,  " 
Sensors  and  Actuators,  vol.  89,  no.  3,  pp.  259-266,  March,  2001. 

26.  Chang  Liu,  "A  magnetic  flexible  membrane  actuator  for  micro  fluid  pumping 
(Invited)  ",  Sixth  International  Symposium  on  Magnetic  Materials,  Processes  and 
Devices,  a  Symposium  at  the  198th  Meeting  of  the  Electrochemical  Society,  Oct 
22-27,  2000,  Phoenix,  Arizona. 


49 


27.  Melvin  Khoo,  C.  Liu,  "Development  of  A  Novel  Micromachined 
Magnetostatic  Membrane  Actuator" ,  The  58th  Annual  IEEE  Device 
Research  Conference,  pp.  109-110,  Denver,  CO,  2000. 

28.  Melvin  Khoo,  C.  Liu,  "A  Novel  Micromachined  Magnetic  Membrane  Microfluid 
Pump",  22nd  Annual  International  Conference  of  the  IEEE  Engineering  in 
Medicine  and  Biology  Society,  2000,  Chicago,  IL. 

29.  BP  Mosier,  JI  Molho  and  JG  Santiago,  Experiments  in  Fluids,  submitted  2001. 

30.  D.  Liu,  M.  Maxey  and  G.E.  Kamiadakis,  A  fast  algorithm  for  particulate 
microflows  in  complex  geometries  Proc.  ASME  Winter  Meeting,  NY,  November 
2001. 

3 1 .  G.  Li  and  N.  R.  Alum,  “A  Lagrangian  approach  for  electrostatic  analysis  of 
deformable  conductors”,  Accepted  for  Publication  in  Journal  of 
Microelectromechanical  Systems,  2001. 

32.  L-H  Lu,  K.  Ryu,  C.  Liu,  "A  Novel  Microstirrer  and  Arrays  for  Microfluidic 
Mixing, "  Fifth  International  Conference  on  Miniaturized  Chemical  and 
Biochemical  Analysis  Systems,  Oct  21-25,  2001. 


50 


APPENDIX  I.  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 


MEMS 

MicroEIectro  Mechanical  Systems 

SPICES 

A  general  purpose  circuit  simulator 

NsKxar 

A  general  purpose  CFD  code  for  simulating  incompressible, 
compressible  and  plasma  flows  in  unsteady  three-dimensional 
geometries 

DSMC 

Direct  Simulation  Monte  Carlo  method  for  fluidic  flows  through 
narrow  channels 

FASTSTOKES 

Precorrected-FFT  accelerated  Stokes  flow  solver. 

CAD 

Computer  Aided  Design 

EOFLOW 

EFT 

Fast  Fourier  Transform 

FASTWAMIT 

Fast  analysis  program  for  offshore  structure  analysis 

MIT 

Massachusetts  Institute  of  Technology 

tspice 

Time  for  membrane  actuation 

tcall 

Time  for  the  next  call 

MEMS 

MicroEIectro  Mechanical  Systems 

V 

Voltage 

Vout 

Output  Voltage 

Rn 

Resistance  value  of  n 

Tn 

Solution  at  n 

Tsso 

Steady  state  solution  of  neural  network 

P 

Pressure 

U 

Flow  velocity 

SEM 

Scanning  electron  Microscope 

Q 

Flowrate 

X 

Heat  Source 

LCP 

Liquid  Crystal  Polymer 

CFD 

Computation  Fluid  Dynamics 

LCP 

Liquid  crystal  polymer 
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